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ABSTRACT 

The  estimation  of  the  frequency  spectrum  of  signals  is  of  widespread 
importance  in  such  diverse  areas  as  seismic  signal  processing,  vibration  mea- 
surements, biomedical  signal  processing,  speech  transmission,  etc.  In  commu- 
nication systems,  the  estimation  of  the  spectrum  of  a signal  is  of  importance 
for  both  demodulation  of  angle  modulated  signals  and  in  frequency  synchroni- 
zation. A commonly  used  technique  for  spectral  estimation  in  stochastic 
signals  is  to  identify  the  parameters  of  an  autoregressive  model  of  the  sig- 
nal and  then  determine  the  spectrum  in  terms  of  these  parameters.  Instead 
of  identifying  the  autoregressive  parameters,  alternative  parameters  may  be 
formulated  which  are  preferable  from  an  identification  viewpoint.  In  this 
investigation,  one  such  set  of  parameters  is  determined  on  the  basis  of  sen- 
sitivity studies.  The  parameters  are  obtained  from  a nonlinear  transforma- 
tion of  the  autoregressive  parameters.  A maximum  likelihood  estimation  scheme 
is  used  to  obtain  appropriate  sequential  algorithms  for  identification  of 
these  parameters.  Several  simulation  examples  are  presented  to  illustrate 
the  applicability  of  the  algorithms  for  spectral  estimation  in  a stationary 
stochastic  process.  In  addition  two  examples  in  communication  systems  are 
presented  to  illustrate  the  applicability  of  the  technique.  The  first  is 
concerned  with  identifying  an  emitter  by  means  of  its  spectral  characteris- 
tics and  involves  adaptively  determining  the  instantaneous  frequency  in  the 
signal.  The  second  example  is  in  FN  code  acquisition  in  a frequency  hopping 
spread  spectrum  system  and  illustrates  the  computational  savings  achieveable 
over  other  comparable  methods  presented  in  the  literature. 
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1. Introduction: 


Interest  in  the  area  of  estimation  of  the  frequency  spectrum  of 
signals  has  grown  significantly  in  recent  years.  Estimation  of  the  fre- 
quency spectrum  has  found  wide  applications  in  seismic  signal  processing, 
biomedical  signal  processing,  speech  transmission,  analysis  of  Doppler 
radar  returns,  vibration  measurements,  etc.  A variety  of  Fourier  trans- 
forms methods  have  been  developed  for  spectral  estimation.  An  alternate 
technique  for  spectrum  estimation  is  to  fit  an  autoregressive  (AR)  or 
autoregressive  moving  average  (ARMA)  model  to  the  process.  For  a signal 
modeled  by  the  AR  process 
N 

y(k)  = l an  y(k-n)  + w(k)  (1) 

n*l 


where  w(*)  is  assumed  to  be  a zero  mean,  Gaussian  white  process  with  vari- 
ance V , the  spectrum  is  given  by 


N 

1-1 

n-1 


a e 
n 


-jwn 


V 

w 


(2) 


so  that  the  problem  of  spectrum  estimation  becomes  one  of  identifying  the 
parameters  of  an  AR  model  of  the  process. 

Several  techniques  are  available  for  the  estimation  of  the  parameters 
in  a general  AR  model  of  the  form  of  Eq.  (1).  The  techniques  basically 
require  the  solution  of  a set  of  equations  involving  the  autocorrelation 
coefficients  and  require  batch-processing  of  the  data.  In  such  techni- 
ques, the  computation  of  the  correlation  coefficients  constitues  a sub- 
stantial part  of  the  total  computational  effort  in  obtaining  the  autore- 
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gressive  parameters. 

In  many  applications,  the  signal  is  not  truly  stationary  and  the 
spectral  content  of  the  signal  may  be  varying.  Typical  examples  are  in 
speech  or  image  processing  or  in  communications.  The  short-time  suectrum 
of  a speech  signal  varies  from  interval  to  interval,  so  that  the  para- 
meters of  the  associated  AR  model  must  be  recomputed  periodically.  In 
image  processing,  the  image  intensity  will  not  be  the  same  over  the  en- 
tire picture.  Thus,  if  an  AR  model  is  used  to  describe  the  image,  the 
parameters  of  this  model  will  not  be  constant  over  the  entire  image.  In 
communication  systems  involving  frequency  or  phase-modulation,  the  wave- 
forms are  typically  characterized  as  having  a narrowband  time-varying 
spectrum.  Demodulation  then  essentially  requires  estimating  the  instant- 
aneous frequency  in  the  spectrum  of  the  received  signal.  Again,  in  spread 
spectrum  systems  involving  frequency  hopping,  the  receiver  is  required  to 
determine  the  instantaneous  frequency  in  the  received  signal.  In  all 
these  cases,  we  therefore  need  to  track  the  coefficients  of  the  AR  model 
adaptively  in  order  to  be  able  to  estimate  the  spectrum  of  the  signal. 

It  is  clear  from  Eq.  (1)  that  once  the  AR  coefficients  {<*n}  have  been 

identified,  the  short-term  spectrum  of  the  signal  y(‘)  can  be  determined 

from  Eq.  (2).  However,  from  the  point  of  view  of  estimating  the  signal 

spectrum,  the  parameters  {a  } of  the  autoregressive  model  may  not  be  the 

n 

best  set  for  identification.  Recent  investigations  have  shown  that 
alternative  coefficients  can  be  associated  with  the  AR  model  which,  for 
a given  data  length,  can  lead  to  better  spectral  estimation.  One  such 
set  of  coefficients  is  the  partial  correlation  coefficients  [l,2j. 

The  present  research  has  been  concerned  with  determining  other  sets 


of  coefficients  which  may  lead  to  better  spectral  estimation  and  to  ob- 
tain algorithms  for  the  identification  of  these  parameters.  Alternative 
coefficients  have  been  determined  which  seem  to  provide  better  estimates 
of  the  short-time  spectrum  using  finite  length  data.  In  the  next  sec- 
tion we  formulate  a set  of  such  coefficients  and  illustrate  the  use  of 
autoregressive  spectral  estimation  techniques  in  two  communications  ex- 
amples. 

2.  Alternative  Coefficients: 

We  now  consider  the  problem  of  formulating  alternative  coefficients 
associated  with  the  model  of  Eq.  (1).  Let  us  denote  by  0 the  (as  yet 
undersignated)  parameter  vector  to  be  identified.  Let  denote  the  set 
of  observations  (y(l),  y(2)  ...  yCk)}.  If  we  characterize  0 as  a non- 
random  parameter,  we  can  use  the  maximum  likelihood  (ML)  technique  for 
identification.  The  ML  estimate  of  0 is  obtained  by  maximizing  the  con- 
ditional density  function  p(Y^  | 0).  We  note  that 

p <Yk I OJ  = P(y(k)|Yk_1,e)p(y(k-l)|Yk_2,e)  ...  p(y (0))  (3) 

so  that  the  log-likelihood  function  is 

k 

£n  p(Y.  |0)  - l In  p (y  (j  ) | Y . .,0)  (4) 

j-1  I - 

Because  of  the  Gaussian  assumptions  on  w(*)»  it  easily  follows  that  y(j) 
is  Gaussian.  From  Eq.  (1),  we  can  write 

N 

Etya ) |Y  0}  — X “.yCj-1) 

~J  i-i 

var{y(j)  | Yj_1,0}  = var{w(j)}  = V w 
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■ 


u 
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The  log-likelihood  function  therefore  is 


k N 

£n  p(Y  |0)  = - | £n(27TVw)  ~ l [y(j)  - l o ylj-t)]2  (5) 

w j=l  i=l 


The  maximum  likelihood  estimate  of  9 is  thus  that  value  of  9 which  mini- 
k N 

mi2es  the  quantity  £ [y(j)  - £ a y(j-i>]  . If  we  define  £(j)  as 

j=l  i-1 


the  residual. 


:(j)  = y(j)  ~ I a y(j-i) 
i=l 


it  is  clear  that  £(•)  represents  the  output  of  a linear  system  with  trans- 
_1  « 

fer  function  H(z  ) = 1 - 2,  ct.  z and  input  y(‘).  The  ML  estimate  9 

i=l  1 

can  then  be  obtained  by  choosing  the  parameters  of  this  system  such  that 
the  square  sum  of  its  output  is  minimized. 

The  preceding  discussion  gives  us  a basis  for  associating  a set  of 
parameters  with  an  AR  process.  We  can  obtain  realizations  of  the  transfer 
function  H(z  h using  different  parameter  sets.  Since  the  identification 
of  the  parameters  involves  minimizing  the  square  sum  of  the  corresponding 
outputs,  we  can  compare  these  parameter  sets  from  an  identification  view- 
point by  considering  the  sensitivity  of  the  output  to  variations  in  the 
parameters.  If  the  sensitivity  is  low,  the  residual  is  not  affected  as 
much  by  variations  in  the  parameter.  This  implies  that  even  if  the  para- 
meters are  not  identified  very  accurately,  H(z  and  hence  $^((0)  can  be 
estimated  well.  Hence  for  a given  data  length  the  parameter  set  with 
lower  sensitivity  should  yield  a better  spectral  estimate. 


One  method  of  comparing  different  realizations  of  H(z  ) is  to  use 
a time-domain  sensitivity  function  defined  as  follows  [3].  Let  "a"  denote 
a parameter  in  the  system  and  s(k,a)  the  corresponding  output  at  the  kth 
stage.  Then  the  sensitivity  function  of  s(.,.)  with  respect  to  "a"  is 
defined  as 


s' 3 00 

a 


3 In  s(k,a) 

3 In  a 

9 s(k,a)  . a 
3 a s(k,a) 


(7) 


Equation  (7)  represents  the  percentage  change  in  the  output  s(k,a)  due 
to  small  perturbations  in  the  values  of  a.  It  is  easy  to  see  from  Eq.  (7) 

g 

that  the  sensitivity  function  S (k)  is  completely  determined  if  the  gra- 

cl 

9 s 

dient  vector  -r—  can  be  computed.  This  can  be  easily  done  by  using  sen- 
da 

sitivity  models  [3,4]. 

Consider  the  linear  system  shown  in  Figure  1 in  which  it  is  desired 
to  find  the  sensitivity  of  the  output  s(k)  with  respect  to  variations  in 
the  parameter  "a".  If  e2(k)  is  added  to  the  output  of  the  gain  element 
"a"  in  the  model,  every  signal  in  the  model  (except  the  output  of  "a") 
is  the  partial  derivative  with  respect  to  "a"  of  the  corresponding  sig- 
nal in  the  system.  This  can  be  seen  by  considering  the  equations  of  the 

system.  From  Figure  1 

k k 

e,(k)  - [ w17(k-j)  u(j)  + l w--(k-j)  e.(j)  + e2(0) 

j-0  j-0 

e3(k)  - ae2(k) 

k k 

e4(k)  = aCk)  - £ w14(k-j)  u(j)  + £ w34(k-j)e3(j)  + e4(0) 

j-0  j-0 


(8) 


5 
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where  represents  the  system  weighting  function  between  the  input  node 
i and  output  node  j with  the  output  of  "a"  open-circuited.  Partial  dif- 
ferentiation of  the  above  equations  with  respect  to  "a"  yields 


3e2(k) 

3a 


3e3(k) 

3a 


3s  (k) 
3a 


k 3e3(j) 

=j=0  W32(k_j) 


3e2(k) 


3a  + e2^ 


k 3e-(j) 

■ I w34(k-j) 

j=0 


3a 


(9) 


Thus  the  gradients  are  obtained  as  the  outputs  from  the  exact  model  with 
zero  input  and  e2(k)  fed  at  the  output  of  the  gain  element  "a". 

When  the  sensitivity  with  respect  to  several  parameters  is  needed, 
it  is  necessary  to  use  as  many  models  as  parameters.  If  the  outputs  of 
these  parameters  end  at  a common  summing  junction  however,  all  the  sen- 
siti  'ity  functions  can  is  determined  simultaneously  by  the  use  of  a single 
model 


3.  Choice  of  Parameter  Set  for  Identification 

We  now  consider  various  realizations  of  the  transfer  function  H(_z 
in  terms  of  different  parameter  sets  and  investigate  their  sensitivity 
properties  with  a view  to  determining  the  best  set  for  identification. 

As  indicated  earlier,  one  set  of  coefficients  is  the  partial  cor- 
relation coefficients.  The  partial  correlation  coefficient  p^,  n = 1 ...  N 
are  obtained  from  the  AR  parameters  by  the  following  set  of  recursions: 

For  n - N,N-1  . . . , 1 


6 


■M 


i 


V 


1 


i = 1,2  ...  n-1 


n-1 


1 - P. 


r n n . 

„ [a.  - p a J 

2 i n n-i 


n-1  n . 

an  = a = 1 
0 o 


n-1 

p = a . 
Hn  n-1 


with 


N 

ai  ~ “i 


PN  ” aN 


i = 1,2  ...  N 


(10) 


T<~  obtain  a realization  in  terms  of  these  parameters,  let  us  define  the 
polynomials 

n 


Fn(z_1)  = l aiZ_i’ 

n i=0  1 


n = s,s-l,  . . . , 1 


(11a) 


and 


G(z-  ) = l 


n -l 
a ,z 
n-i 


(lib) 


i=0 


Now  Eq . (10)  can  be  rewritten  as 


n n-1  , n-1 

a,  = a.  + p a 
i i n n-i 


i = 1 n-1 


(12) 


It  then  follows  by  direct  substitution  of  Eq.  (12)  in  Eq.  (11a)  and 
Eq.  (lib)  that 


F (z-1)  = F (z-1)  + p G (z_1) 
n n-1  n n-1 


(13) 


and 


’Vrl>  ■ Vi'2'1)  + '.vf1* 


(14) 


For  any  input  X(z)  therefore,  we  can  determine  the  output  of  the  system 

with  transfer  function  F (z  F)  in  terms  of  the  outputs  F , (z  '*’)X(z)  and 

n n— 1 


7 


G . (z  ^)X(z)  by  use  of  Eqs . (13  and  (14).  Recursive  use  of  these  equations 
n-i  1 

then  yields  the  lattice  realization  of  Fig.  2 for  the  all  zero  transfer 
function  H(z  ^) . The  initial  values  for  starting  the  recursion  are 

Fq (z  1)  = 1 and  GQ(z  1)  - z 1 


The  coefficients  p have  also  been  referred  to  as  reflection  coefficients, 
n 

The  name  comes  from  transmission  line  theory  in  which  p^  is  considered 

as  the  reflection  coefficient  at  the  boundary  between  two  sections  of 

transmission  line  with  characteristic  impedances  Z and  Z , . . p is  then 

r n n+1  n 

given  by 


Z , . - Z 
n+1  n 

~ Z . + Z 
n+1  n 


1 < n < N 


(15) 


The  transfer  function  H(z  1)  can  then  be  considered  as  that  of  a sequence 
of  these  sections  with  impedance  ratios  given  by 


1 * Pn 
1 - Pn 


1 < n < N 


(16) 


The  same  phenomenon  holds  for  any  type  of  situation  where  there  is  a 
plane  wave  transmission  with  normal  incidence  on  sections  with  different 
impedances.  For  example,  for  the  case  of  an  acoustic  tube  with  n sec- 
tions of  varying  cross-sections,  the  Impedance  ratios  reduce  to  the  in- 
verse ratio  of  the  consecutive  cross-sectional  areas.  An  example  is  the 
acoustical  tube  model  of  the  vocal  tract  used  in  speech  processing  appli- 
cations [5].  We  can  use  as  a second  set  of  parameters  to  be  used  in  spec- 
tral estimation,  the  set  (S  } defined  as 

n 


S 

n 


1 - Pn 
1 + Pn 


n - 1,2  ...  N 


(17) 


we 


In  order  to  obtain  a realization  of  H(z  1)  in  terms  of  {S^} , 
substitute  Eq.  (17)  into  Eqs.  (13)  and  (14)  to  obtain 

<l  + V V2"1' ' « + s„>  W2"l>  + n - V W2_1> 

(18) 

(1  + Sn)  Gn(2_1)  = z-1  [ (1  + Sn)  Gn_1(z“1)  + (1  - SQ) 

Fn_i(z"1)]  (19) 

One  stage  of  the  realization  is  shown  in  Figure  3. 

In  order  to  compare  the  two  parameter  sets  with  regard  to  their  sen- 
sitivity properties,  we  use  sensitivity  models  as  discussed  in  the  pre- 
vious section.  The  sensitivity  model  for  the  parameters  (pn)  is  shown  in 
Figure  4.  A similar  model  can  be  obtained  for  the  parameters  {Sn)  also  [6]. 

In  order  to  compare  the  two  structures  to  determine  which  is  less 
sensitive,  we  use  as  a measure  of  sensitivity,  the  following: 


k=kf 

l 

k=0 


N 


l 

i=l 


(20) 


It  is  easy  to  verify  that  the  time-domain  sensitivity  measure  defined  by 
Eq.  (20)  is  equivalent  to  the  frequency  domain  sensitivity  measure  de- 
fined in  [7 J , as 


3s  _ Lim 
3pi  Ak^-K) 


ir 


I log  $y(di,w) 


-IT 


- log  $y  (d_^+Ad^  ,oj)  | doj 


(21) 
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_ 


t 


I 


I 

I 

I 

\ 

where  (•  ,uj)  is  the  spectrum  of  the  signal  y(*)  [see  Eq . (2)], 
and  the  quantity  in  brackets  is  the  spectral  deviation  due  to  a per- 
turbation in  the  i^  coefficient. 

| 

Several  examples  were  run  to  compare  the  sensitivities  of  tnc  two 
realizations.  We  present  here  a typical  example  in  which  H(z  was 
chosen  as 

H(z_1)  = 1 - 2.8775  z'1  + 2.5673  z~2  - *7897  z-3 

The  corresponding  parameters  obtained  from  Eqs.  (10)  and  (17)  are 

px  = 0.9657,  p2  = 0.7820,  P;J  = 0.2683 

S1  = 0.0175,  S2  = 0.1223,  S3  = 0.5825 
Figures  5 and  6 show  the  sensitivities  with  respect  to  the  two  sets  of 
parameters  when  the  input  to  the  sensitivity  model  is  a unit  step  in- 
put. As  can  be  seen,  the  sensitivities  with  respect  to  the  parameters 

{Sn>  are  less  than  with  respect  to  the  {pnL  This  is  also  illustrated  in 
Table  I which  shows  the  sensitivity  measure  of  Eq.  (20)  for  the  example. 

In  general,  the  examples  show  that  when  the  zeros  of  H(z  ^)  are  close  to 
the  unit  circle,  the  sensitivity  with  respect  to  {Sn}  is  better  than 
with  respect  to  {pn}. 

While  the  discussions  above  show  that  a realization  of  H(z  3)  using 

the  {S  } parameters  does  indeed  exhibit  lower  sensitivity,  these  coeffi- 
n 

cients  are  not  really  suitable  from  an  identification  viewpoint  since 

some  of  the  nice  symmetry  properties  associated  with  the  {p  } parameters 

n 

no  longer  hold.  For  example,  p^  in  the  range  (-1,0)  yields  a value  for 
in  (l,00),  whereas  p^  in  (0,1)  corresponds  to  Sr  in  (0,1),  so  that  the 
range  of  is  (0,°°).  However,  if  we  define  the  coefficients  as 

10 
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Ln  - logSn  n = 1,  ...  N (22) 

it  is  clear  that  we  will  retain  the  symmetry  properties  of  {pp}  while 
r obtaining  lower  sensitivity.  For  these  reasons,  the  parameter  vector 

0 was  chosen  to  be 

2 - [ 4 l2  • • • h ]T 

where  N defines  the  order  of  the  model. 

4.  Identification  Algorithm 

i We  now  consider  the  estimation  of  the  parameters  L^.  In  order  to 

/\ 

obtain  the  maximum  likelihood  estimates  L^,  we  follow  the  procedure  in  [4] 
and  minimize  the  output  £(•)  of  the  system  with  transfer-function 
1 N 

H(z  ; = 1 - l a z . We  use  a realization  of  this  transfer  function 
! i-1 

in  terms  of  the  parameters  (L  } obtained  by  substituting  Eq.  (22)  into 

n 

A 

Eqs.  (18)  and  (19)  and  seek  an  estimate  {L^}  such  that  the  square  sum  of 
the  output  is  minimized.  The  realization  is  shown  in  Figure  7.  From  the 
figure,  the  equations  describing  the  realization  c’an  be  written  down  as 


n+l,t  . 

ef  (k) 


_n„  „ , l-eLn+l  _n„  , x 
e,(k)  + = e,  (k-1) 

f l+eLn+l  b 


1,2, 


(23) 


n+1 

S 


00 


£j(k-D  + 


l-eLn+l 

l+e^n+l 


ej(k) 


(24) 


with  £°(k)  - *;°(k)  “ y(k) . The  estimate  obtained  by  minimizing  the 

square  sum  of  the  residuals  {tjn+^(*)}  are  given  by 


11 


Ww 


£n 


I []ef(i)S(1'1)+{Gb(1-1)}^] 

j1[{^(i“1>}2'e”(i)^(i“l)j 


(25) 


where  * refers  to  estimated  values. 
Minimizing  {e^+^(*)}  similarly  yields 


Ln+l(k)  “ 


l I eJUJejJCi-D+teJU)}' 


Z [tej 

i-1  L 


(i)}2-e"(i)e£(i-l) 


(26) 


Since  the  output  y(k)  is  assumed  to  be  stationary,  it  can  be  shown  by 

k _ k 

direct  computation  that  for  large  values  of  k,  £ {E"(i-l)}Z  = £ {e£(i)}Z 


i-1 


i-1 


We  can  then  combine  Eqs.  (25)  and  (26)  to  get 


Wk)  ” ln 


l {E*(i)+e£(i-l)}2 


i-1 


I {^(D-ePU-I)}' 

i-1  D 


(27) 


Equation  (26)  can  be  put  into  sequential  form  to  yield 


Wk>  “ W1^  + 


1 + 


£n 


1 + 


{e"(k)-re£(k-l)}2 

J^{e"(i)+e"(i-l)}2 


{e^(k)-e^(k-l)}2 

j (ej(i)-ej(i-l)}3 


(28) 


Equations  (23),  (24)  and  (28)  constitute  a complete  set  of  algorithms  for 


identifying  the  parameters  {L  }. 

n 
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5 . Simulation  Examples 


Several  examples  were  simulated  to  illustrate  the  use  of  the  identi- 
fication algorithm.  These  results  are  sunmarized  in  Table  II.  Example  1 
is  an  AR  process  with  real  poles  at  0.9  and  -0.9.  The  estimated  parameters 
are  within  one  percent.  Example  2 is  an  AR  process  with  a pair  of  complex 
poles  at  a radius  of  0.95  and  the  estimated  parameters  are  within  2.32  per- 
cent. Examples  3 and  4 are  AR  processes  with  4 real  poles.  Note  that  the 
percentage  error  in  the  last  estimated  parameter  becomes  larger  as  the 
poles  approach  the  unit  circle.  Examples  5 to  7 are  AR  processes  with 
complex  poles.  As  before,  the  errors  tend  to  increase  as  the  poles  ap- 
proach the  unit  circle.  In  every  case  a data  length  of  5000  samples  was 
used.  The  initial  estimates  of  the  parameter  were  chosen  as  zero. 

The  estimate  of  the  spectrum  $^((0)  can  be  obtained  from  Eq  (2).  In 
order  to  directly  relate  the  spectrum  to  the  parameters  {L^},  we  can  use  a 
realization  of  the  transfer  function 


G(w) 


N 

1-  7 a e 

i n 

n=l 


-jeon 


(29) 


in  terms  of  these  parameters,  as  shown  in  Fig.  8 [9].  This  realization  is 
obtained  from  Fig.  7 by  reversing  the  direction  of  signal  flow.  A typical 
example  is  shown  in  Fig.  9 which  plots  the  true  and  estimated  spectra  for 
Example  6 in  Table  II.  As  can  be  seen  from  the  figure,  the  estimate  follows 
the  true  spectrum  closely. 

In  addition,  two  typical  examples  of  spectral  estimation  in  communica- 
tion systems  were  investigated.  The  first  example  involved  identifying 
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emitters  by  their  spectral  characteristics.  The  problem  is  to  determine 
the  instantaneous  frequency  of  the  largest  amplitude  sinusoid  within  the 
portion  of  the  spectrum  of  interest  [10],  The  model  is  assumed  to  consist 
of  a single  large-amplitude  sinusoid  (the  signal)  plus  a number  of  smaller 
amplitude  sinusoid  (interference).  The  observations  were  assumed  to  be 
corrupted  by  additive  noise  n(k). 

The  data  was  generated  using  the  following  model: 

x(k)  = 1000  cos(61(k))  + 100  cos(02(k))  + n(k),  k = 0,1,2,  ....  511 

where 

ex(k)  = 340  it  kAt  - 24  cosCj1  kAt  - |)  + 2 

02(k)  = 320  tt  kAt  - 40  cos(2t r kAt)  -2 

At  = 0.00125 

It  is  desired  to  determine  the  instantaneous  frequency  every  0.04  sec,  where 
L the  instantaneous  frequency  is  defined  as  ~ and  $ is  the  argument 

of  the  sinusoid.  The  noise  n(k)  is  a sequence  of  pseudo-unif ormly  distri- 
buted random  numbers  in  the  range  -100  n(k)  100.  The  spectrum  of 

interest  is  in  the  range  from  130  to  270  Hertz.  The  sampling  rate  was 
chosen  to  be  800  samples/sec. 

The  power  spectrum  was  estimated  for  every  32  samples.  The  estima- 
tion algorithm  was  used  repeatedly  on  the  same  set  of  data  points  till 

convergence  was  obtained.  At  each  iteration  the  last  estimated  value 
from  the  previous  estimate  was  used  as  the  initial  estimate.  It  was  found 
that  while  the  algorithms  had  to  be  iteratively  applied  30  times  over  the 


I 

w 

f 

1 

I 

I 
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first  data  set  of  32  samples,  only  4 iterations  were  required  for  subse- 
quent intervals  to  achieve  reasonable  convergence  of  the  estimates.  The 
results  are  tabulated  in  Table  III  for  model  orders  2,4  and  6. 

The  second  communication  example  was  taken  from  spread  spectrum 
systems.  Spread  spectrum  systems  have  found  increasing  applicability  in 
areas  such  as  remote  piloted  vehicles,  anti-jam  systems,  packet  radio  net- 
work, multiple-access  systems,  etc.  A spread  spectrum  signal  is  generated 
by  modulating  a data  signal  onto  a wide  band  carrier  so  that  the  resultant 
transmitted  signal  has  bandwidth  much  larger  than  the  data  signal.  In  \ 

the  frequency  hopping  (FH)  spread  spectrum  system  an  n-stage  PN  sequence 
generator  is  used  to  determine  the  carrier  frequency  over  each  signalling 
interval  (called  a chip  interval) . The  carrier  is  then  modulated  by  the 
data  before  transmission.  At  the  receiver,  the  signal  is  demodulated  by 
correlating  it  with  a signal  generated  locally.  For  proper  demodulation, 
the  frequency  of  the  local  signal  must  be  the  same  as  that  of  the  carrier. 

It  follows  that  a primary  requirement  in  such  receivers  is  to  set  the  ini- 
tial state  and  maintain  the  phase  of  the  locally  generated  signal  in  syn- 
chronism with  the  received  PN  signal.  In  this  example  we  consider  the 
application  of  the  spectral  estimation  technique  discussed  earlier  to  the 
problem  of  PN  code  acquisition  at  the  receiver  of  a frequency  hopping  sy- 
stem [ ll] . For  purposes  of  clarity  we  will  assume  that  the  data  consists 
of  a single  bit  over  each  chip  interval. 

A simple  block  diagram  of  an  FH  transmitter  is  shown  in  Figure  10.  An 
n-stage  linear  feedback  shift  register  generates  a PN  sequence  of  length 
2n-l.  The  output  from  the  shift  register  stages  is  used  as  the  binary  in- 
put bnbn_^. . *^2^1  to  a frequency  synthesizer  which  generates  a signal  at 
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the  frequency  corresponding  to  the  binary  number  b^b^^. . .b^b^,  where  b^ 
is  the  data  bit.  The  frequency  hopping  rate  can  be  either  equal  to  or  a 
multiple  of  the  data  rate.  We  will  assume  here  that  the  PN  code  bit  rate 
is  equal  to  the  data  bit  rate.  The  synthesizer  output  is  suitably  ampli- 
fied and  upconverted  to  the  required  frequency  for  transmission.  The 
transmitted  signal  can  be  expressed  as 

s(t)  = J"2P  £°°  cos(27Tfkt  + <|>k)  rect 

lc=-oo 

where  f^  is  the  hopped  frequency  corresponding  to  the  (n+l)-bit  binary 
number  during  the  kth  chip  interval  where  each  interval  is  of  duration  T. 
The  phase  4>k  is  assumed  to  be  uniformly  distributed  between  0 and  2ir  for 
non-coherent  frequency  synthesizers.  The  instantaneous  spectrum  of  the 
transmitted  signal  thus  is  a single  discrete  line. 

The  operation  of  the  FH  receiver  can  be  easily  explained  with  the 
aid  of  the  block  diagram  in  Figure  11.  A shift  register  is  used  to  gener- 
ate a replica  of  the  PN  sequence  used  at  the  transmitter.  The  n-bit  binary 
number  corresponding  to  the  shift  register  state  is  used  as  input  to  the 
frequency  synthesizer  which  generates  one  of  the  possible  2n  frequency 
signals.  The  received  signal  is  appropriately  down-converted  and  multi- 
plied with  the  locally  generated  signal  resulting  in  an  intermediate  fre- 
quency signal  that  is  frequency  shift  keyed  by  the  data,  so  that 

00 

r'Ct)  ■ A l cosC2irfIF  + 2irAfk  + 1^)  rect(t~k— ) + n'  (t) 
k=-°° 

where  Afk  corresponds  to  the  kth  data  bit,  and  n'(t)  represents  the  noise 
term.  This  signal  is  then  demodulated  by  conventional  techniques. 
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Before  a meaningful  FSK  demodulation  can  be  attempted,  the  locally 
generated  PN  sequence  must  be  synchronized  to  the  received  sequence. 

This  synchronization  is  a two-step  process:  initial  acquisition  and 
tracking . 

It  is  clear  from  the  previous  example,  that  the  identification  al- 
gorithm presented  in  this  report  can  be  used  to  identify  the  instantaneous 
frequency  in  the  spread  signal  and  thus  provide  tracking.  However,  ini- 
tial acquisition  can  be  time  consuming  if  a serial  search  is  conducted 
over  all  possible  phases  of  a long  PN  sequence.  The  problem  of  PN  code 
acquisition  at  the  receiver  of  a FH  system  is  aimed  at  rapid  acquisi- 
tion of  the  initial  state  of  the  transmitter  PN  code  generator.  This  ini- 
tial state  is  used  in  the  receiver  PN  code  generator  which  is  a hardware 
replica  of  the  transmitter  PN  code  generator.  After  PN  code  acquisition 
the  receiver  PN  code  generator  is  locked  on  to  the  FH  receiver.  One 
such  PN  code  acquisition  scheme  ASEAT,  Autoregressive  Spectral  Estima- 
tion Acquisition  Technique  is  considered  in  [ll].  The  technique  essen- 
tially consists  of  estimating  n successive  transmitted  frequencies  and 
hence  n successive  PN  code  generator  states.  The  most  significant  bits 
(MSB)  of  these  n successive  PN  code  generator  states  are  shifted  into 
another  shift  register  which  at  the  end  of  this  n stage  operation  will 
contain  the  initial  PN  code  generator  state  to  a high  degree  of  accuracy. 
In  effect  this  assumes  the  structure  of  the  PN  code  generator  to  be  a 
shift  register  with  a feedback  at  the  least  significant  bit.  Thus  the 
acquisition  logic  uses  the  MSB  of  the  estimated  PN  state  at  each  stage 
and  hence  acquisition  time  is  of  the  order  of  nT  seconds.  If  the  es- 
timated state  is  accurate  to  k most  significant  bits  then  the  acquisition 


17 


> 


logic  can  be  modified  to  given  an  acquisition  time  of  the  order  of 
(n-k+l)T  seconds. 

The  spectral  estimation  algorithm  presented  in  the  earlier  sections 
of  this  report  was  applied  to  the  problem  of  PN  code  acquisition  in  the 
FH  system  described  above.  For  purposes  of  comparison  with  the  results 
in  [ll],  the  system  was  simulated  using  a data  rate  of  10  Kbps  and  with 
n = 9.  The  number  of  hopping  frequencies  used  was  equal  to  1022.  The 
9-stage  PN-code  generator  is  a maximal  length  sequence  generator,  gene- 
rated by  the  primitive  polynomial 

4 9 

P(x)  = 1 + x + x 

The  PN  code  generator  is  realized  using  a shift  register  and  an  exclusive 
or  gate  as  shown  in  Fig.  11.  The  initial  and  successive  states  of  the 
PN  code  generator  at  the  transmitter  are  406,300,89,179,359,207,414,316, 
120,241,  ...  etc.  The  data  bit  is  augumented  as  the  least  significant 
bit  to  the  9-bit  PN  generator  state  and  the  resulting  number  is  used  to 
generate  the  transmitted  signal  with  one  of  the  1022  uniformly  spaced  fre- 
quencies lying  in  the  interval  O.Itt  < < 0.78tt.  The  received  signal  is 

obtained  by  adding  zero-mean  white  gaussian  noise  such  that  the  carrier 
to  noise  ratio  is  -12db.  The  estimation  algorithm  is  used  to  estimate 
the  transmitted  frequency  and  hence  the  state  of  the  PN-generator . 

The  successive  PN  generator  states  and  its  estimates  are  summarized 
in  Table  IV.  The  estimates  for  model  order  4 are  comparable  to  the  results 
in  [ll]  which  uses  a model  order  16  and  hence  a large  reduction  in  model 
order  is  obtained.  Since  the  acquisition  logic  in  [ll]  uses  only  the 
MSB  of  the  estimated  state  in  each  successive  state,  from  Table  V it  can 
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be  seen  that  even  a model  order  of  2 can  be  used  without  any  acquisi- 
tion error,  thus  reducing  the  computational  efforts  considerably. 

As  suggested  earlier  the  accuracy  of  the  estimation  scheme  allows 
the  use  of  more  bits  of  each  estimated  state.  Successive  estimates  of 
a particular  bit  of  the  initial  state  can  be  passed  through  a majority 
decoder  to  further  reduce  the  error  probability  in  the  estimation  of 
the  initial  state.  An  acquisition  scheme  using  an  estimator  of  model 
order  4,  6 most  significant  bits  of  the  estimated  state  at  each  state 

and  a majority  3 decoder  for  each  bit  of  the  initial  state  to  be  esti-  ' 

mated  is  shown  in  Table  V . The  acquisition  logic  acquires  the  initial 

state  to  be  406  which  in  fact  was  the  transmitted  initial  state.  The 

lock-in  time  for  the  acquisition  logic  is  4T  seconds  as  compared  to  9T 

seconds  for  the  acquisition  logic  in  [ll]. 

6.  Conclusions 

This  report  has  considered  the  estimation  of  parameters  in  an  auto- 
regressive model  of  a stationary  stochastic  process.  A set  of  parameters 
which  are  more  suitable  from  an  identification  point  of  view  were  deter- 
mined on  the  basis  of  sensitivity  studies.  The  parameters  are  obtained 
from  a nonlinear  transformation  of  the  autoregressive  parameters.  A 
maximum  likelihood  estimation  scheme  was  used  to  obtain  appropriate  sequ- 
ential algorithms  for  identification  of  these  parameters.  Several  simu- 
lation examples  are  presented  to  illustrate  the  applicability  of  the  al- 

gorithms  for  spectral  estimation  in  a stationary  stochastic  process.  | 

One  of  the  original  purposes  of  the  investigation  was  to  determine 
the  applicability  of  the  spectral  estimation  technique  to  the  problem  of 
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demodulation  of  digital  modulated  signals.  With  this  in  view  two  examples 
in  communication  systems  were  considered.  The  first  was  concerned  with 
identifying  an  emitter  by  means  of  its  spectral  characteristics  and  in- 
volved  adaptively  determining  the  instantaneous  frequency  in  the  signal. 

The  second  example  was  in  PN  code  acquisition  in  a frequency  hopping  spread 
spectrum  system.  The  simulation  results  show  the  applicability  of  algorithms 
to  this  problem. 
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TABLE  II 


Order  of 


Example 

No. 

the  AR 

process 

Actual  Locatloh 
of  Poles 

Actual  {L . } 

i 1 

Estimated 
after  5000 

1 

2 

-0.9 

0.0 

-0 . 058299A 

2.25A058 

2.2753185 

2 

2 

0. 95e^Z 

1.759035A 

1.75A8A3 

-2.9710717 

-2.9023A7A 

3 

A 

0.5 

A . 1A0AA99 

A . 11A2841 

0.55 

-2.6003106 

-2.6100783 

0.6 

1.1321A68 

1.0957899 

0.65 

-0.2153282 

-0. 166A129 

A 

A 

0.8 

7 . A872535 

7 . 57269A3 

0.85 

-5 . 78615A6 

-5.7197308 

0.9 

3.830277A 

3. A6A6811 

0.95 

-1.32915 

-0.A76785 

i-jiL 

5 

A 

0 . 7e  3 6 

2.660932A 

2.669A032 

±JL 

-2.06076A6 

-2.0663978 

0 . 5e  3 A 

0.8972335 

0.8581285 

-0 . 2A62367 

-0.2A15095 

+ .7T 

6 

A 

0 . 8e  3 6 

2 . 3'107328 

2 . 3AA77A7 

+ .TT 

-2. id5060A 

-2. 283A3A8 

0 . 7e_33 

1.165021 

1.1353965 

-0.6A90662 

-0.63825A5 

7 

A 

0 . 95e  3A 

1 . 82279A1 

1. 822A967 

41  trr 

-A. 1579523 

-A. 0788953 

0 . 7e"^A 

1.5375355 

1.528A85 

-0.9A99867 

-0.8927A5 
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FIG.  2.  REALIZATION  OF  H(z~ 
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FIG.  6:  SENSITIVITY  FUNCTIONS  OF  {S  } 
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ACTUAL  AND  ESTIMATED  SPECTRA  FOR  AN  EXAMPLE 
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FIG.  10.  BLOCK  DIAGRAM  OF  AN  FH  TRANSMITTER 
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FIG.  11.  BLOCK  DIAGRAM  OF  AN  FH  RECEIVER 


v 


